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Abstract. The modeling of the elastic properties of disordered or nanoscale solids requires the foundations of the 
theory of elasticity to be revisited, as one explores scales at which this theory may no longer hold. The only cases for 
which microscopically based derivations of elasticity are documented are (nearly) uniformly strained lattices. A micro- 
scopic approach to elasticity is proposed. As a first step, microscopically exact expressions for the displacement, strain 
and stress fields are derived. Conditions under which linear elastic constitutive relations hold are studied theoretically 
and numerically. It turns out that standard continuum elasticity is not self-evident, and applies only above certain spa- 
tial scales, which depend on details of the considered system and boundary conditions. Possible relevance to granular 
materials is briefly discussed. 



PACS. 46.25. Cc Static elasticity: theoretical studies - 61.43. -j Disordered solids - 
of nanoscale materials - 83.80.Fg Granular solids 
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1 Introduction 

It is quite surprising that the existing microscopic justification 
of the time-honored theory of elasticity, which has been thor- 
oughly researched in a variety of disciplines, is limited to lat- 
tice atomic configurations [Q], Classical continuum elasticity 
theory has been applied to a large variety of systems, including 
granular materials [|[. In recent years the same theory has been 
applied for the description of elastic properties of micro- and 
nano-scale systems (e.g., [j|]). It is a-priori unclear whether this 
theory applies at such small scales. 

The study presented below shows that the justification of 
elastic theory based on a microscopic picture is not entirely 
straightforward. As expected, one finds that linear continuum 
elasticity is valid on sufficiently large scales. Another result 
is that like in granular matter, one observes force chains in 
strained elastic systems (also observed in Qj). These chains 
are not "visible" in the corresponding stress field. Classical me- 
chanics is assumed throughout this paper. The case of isostatic 
systems, which has received considerable interest in the litera- 
ture, is not specifically addressed here. 



2 Coarse Graining and Constitutive 
Relations 

2.1 Preliminaries 

Consider a system of particles (indexed by {{}) whose masses, 
center of mass positions and velocities at time t are given by 



{m,; r,- (t) ; v,- (t) }. Following [||] define the coarse grained mass 
density at position r and time t as 



i 

Similarly, define the momentum density as 



(1) 



(2) 



where </>(R) is a normalized non-negative coarse graining func- 
tion (with a single maximum at R = 0) of width to, the coarse 
graining scale. Unlike in [§], here only spatial, and not tempo- 
ral coarse graining, is invoked. Upon taking the time derivative 



of Eqs. (|1| 
ulations 



and performing straightforward algebraic manip- 
one obtains two of the equations of continuum me- 



chanics. Eq. (|1|) yields the equation of continuity: 

p(r, t) = - V ■ p(r, t) = - V [p(r, t)V(r, <)] , (3) 

where p = and the coarse grained velocity field is defined 
by V(r, t) = p(r, t)/p(r, t). From Eq. (Q) one obtains the mo- 
mentum conservation equation: 



Pa(r,t) 



_d_ 

dr 







[p(r, t)V a (r, t)V p (r, t) - a a p (r, t)} , (4) 



where Greek indices denote Cartesian coordinates. 

Define fij a (t) to be the a-th component of the force ex- 
erted on particle i by particle j {j ^ i) at time t (assuming pair- 
wise interactions), = — r 3 , and the fluctuating velocity 
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of particle i: v^(r, t) = v,(t) — V(r, t). With these definitions, 
the stress tensor, a a p, is given by the following expression: 



<T a p(r,t) =-- 



E fi0Oi{t)njfi{t) (5) 



x / dscf>[r - Ti(t) + STij(t)] 
Jo 

■ m i v 'ia( r > t)v' i0 (r, t)4>[r - Ti(t)], 



The first term in Eq. (||) is commonly referred to as the "con- 
tact stress" (or "collisional stress"), while the second term is a 
kinetic contribution (the "kinetic stress" or "streaming stress"), 
which vanishes for quasi-static deformations. 

The energy conservation equation can be obtained in a sim- 
ilar way. Assume, for simplicity, that the forces are derived 
from a potential function: fy = —'Vi'P (r.y), with obvious no- 
tation. Define the energy density as: 

i 

+ \ E *(r«(t))^[r-ri(t)]. (6) 

Application of a time derivative to e(r, t) and a rearrangement 
of terms yields the energy conservation equation: 

d 

^( r ' *) = ~7T~ ^( r ' *) e ( r ' *) ~ VaCTapir, t) + cp(r, t)} , 

(7) 

(8) 

«^( r »*M r - r *(*)] 



where the heat flux, c, is given by: 



-I- t E K ( r > *) + u i« ( r > *)] fii« (t)n jf ) (t) 



x / ds0[r - r;(i) + sr^i)] 
Jo 



2.2 Displacement and Strain 

Following elementary continuum mechanics, consider a mate- 
rial particle whose initial (Lagrangian, at time t = 0) coor- 
dinate is R. Its (Eulerian) coordinate at time t is denoted by 
r(R, t). The corresponding (Lagrangian) displacement field is 
given by u La (R, t) = r(R, t) — R. The material particle's ve- 
locity is V La (R, t) = <9u La (R, t) /dt. It therefore follows that 
u La (R. ,t) — J * V La (R, t')dt' . Using the definitions presented 
in Sec. ETL one obtains: 



u La (R,i)^ /' 
Jo 



D i m < v i (tO^[r(R,t / )-r i (tQ] 
£,m^[r(R,i')-r^')] 



di'. (9) 



The macroscopic displacement field u is history dependent, 
i.e., the displacement at time t depends on the trajectories of 
the particles from t = to t. However, noting that = Vi, 
where = rj(i) — rj(0) is the displacement of particle i, and 
invoking integration by parts in Eq. (fy, one obtains (in the Eu- 
lerian representation): 



u(r, t) 



u lin (r, t) 



(10) 



x^[r-r 4 (t')] 



Y,m l v' lfj (r 1 t')u' la (r 1 t') 



where 



u lm (r,t) 



J2imiUi(t)^[r - Tj(t)] 
E,- TO i^[ r - r i(*)] 



(11) 



(12) 



and u<(r, i) = u % (t) - u lin (r, i). 

It is claimed that u lln represents the displacement field rel- 
evant to linear elasticity, i.e., the error is nonlinear in the strain. 
To this end, let e* = max,, ^"'~ u ^ where j represents the near- 
est neighbors of i and ay is the distance between the particles 
i and j. Let e* = max^ e*. Linear elasticity is a theory which is 
linear in e*. Following the above definitions, 



Ej m i i u ia(t) - %<*(*)] ~ "",(*)] 

J2 k m k (f>[r -r k (t)} 



(13) 



hence \u' ia \ < C±e*, where C\ depends on the coarse grain- 
ing scale, and \v' ic \ < C2£*/r, where r represents the typ- 
ical time scale on which the {u^} change (in the quasistatic 
limit,r — > oo, while e* remains finite). The integrand in Eq. jlo[ ) 



is thus O 



. As a matter of fact, it is easy to show that the in- 



tegral on the right hand side of Eq. ( |10[ ) is bounded from above 
by O (e* 2 ). Note that ^ v' it since is defined with re- 
spect to the linear displacement u hn and not the exact displace- 
ment, u. A useful feature of the linear displacement, u hn (r, t), 
is that, unlike the exact displacement field [Eq. (|To|)], it depends 
only on the microscopic displacements at time t. 

The linear strain field is e a/ g(r, t) = -| ^"^ ^ + ^§r 
Notice that e a p = O (e*). Following the above arguments, this 
expression can be replaced by: 



dr a 



(14) 



the error being O (e* 2 ). 

It is interesting to compare the above results with some pre- 
viously defined heuristic calculations of the strain field. The 
mean field strain (e.g. [||]) is based on the assumption that 
the relative particle displacements are described by the macro- 
scopic strain, i.e., iiy Q (r,i) = eapr^p where = u.; — Uj 
(an affine deformation). An improvement of this method, which 
enables a local evaluation of the strain field, is provided by the 
"best fit" hypothesis [Q], whereby the rms difference between 
the actual relative displacements and the above mean field ex- 
pression for them is minimized in a given volume, to produce 
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a 'best' strain field. The mean field approaches are in general 
inconsistent with local force equilibrium (except for homoge- 
neous deformations of lattice configurations). Therefore, the 
mean field (or best fit) strain constitutes an uncontrolled ap- 
proximation of the strain field. The difference between the ex- 
act strain field, the linear strain as given by Eg. (|l4|), and the 
best fit approximation are demonstrated in Fig. |]. In this figure, 
the above three fields are presented for a one dimensional chain 
of 1000 point particles connected by linear springs of random 
rest lengths (mean rest length a = 1, relative standard deviation 
0.29) and random spring constants (relative standard deviation 
0.13), the global strain of the system being e = 0.05. To this 
end, a Gaussian coarse graining function cf>(r) = j^e"" r "'"' 
with w = 50a has been employed. In the calculation of the best 
fit strain, the fluctuations with respect to the mean field are 
weighted using the same coarse graining function as in the ex- 
act formulation. As seen in Fig. n , the linear strain is very close 
to the exact strain, whereas the best fit provides quite a poor 
approximation (for the case of equal spring constants, the best 
fit method yields a uniform strain, while the correct strain is 
space dependent H). In general, heuristic approximations for 
the strain field may result in inaccurate constitutive relations. 



0.052 
0.051 



0.05 

e(x) 

0.049 
0.048 
0.047 
0.046 





Exact 




Linear Strain 




- - - Best Fit 


\ 




\ - - - ' .X^- 
— \ y ^^v* 













400 



450 



500 550 
X 



600 



650 



Fig. 1. The Eulerian strain e(x) vs. position x in the central region of a 
linear chain of particles connected by linear springs of random spring 
constants and rest lengths, calculated by three methods (see text). 



2.3 Stress-Strain Relation 

It is not a-priori clear that the stress field can be expressed as 
a linear functional of the strain field, even for small deforma- 
tions. Each of these two macroscopic fields, cf. Eqs. (p|Jl4]), is 
a different average of microscopic entities. Once averaging is 
invoked to produce macroscopic fields, on cannot deduce the 
microscopic entities from these fields. This is one of the funda- 
mental flaws of the mean field approaches, which rely on such a 
deduction. Below, an exact method for obtaining linear elastic- 
ity, which demonstrates the above problem and highlights the 
scale limitations of this theory, is outlined. 



Consider, for sake of simplicity, only systems with pair- 
wise interactions. In order to develop linear elasticity one can 
assume, without loss of generality, harmonic interactions: 



*(r ij {t)) = ±K ij QT ij \-l i j) S 



(15) 



where kj is the equilibrium separation of particles i and j. The 
force on particle i exerted by particle j is given, to linear order 
in the relative particle displacements, uy, by: 



-K, 



(16) 



where the superscript denotes the reference configuration, 
in which all particle pairs are at their equilibrium separation 
(|r°- | = lij), i.e., an unstressed configuration (prestressed states 
are not considered here). 

It is apparent that even in this case, the microscopic ex- 
pressions for the contact stress [the first term in Eq. (Bh] and 
the strain [Eq. (p^j)] are not manifestly proportional. Therefore 
macroscopic elasticity is not a-priori obvious. To see how it still 
comes about, substitute Eq. (|l6|) in Eq. (|j). The contact stress, 
to linear order in {u,j }), is: 



ij 

x f d S 0[r-r° + S r°.] . 
Jo 



Consider a volume f2, whose linear dimension, W, is much 
larger than the coarse graining scale, w, and let r be an in- 
terior point of fl which is 'far' from its boundary. Let upper 
case Latin indices denote the particles in the exterior of fl 
which interact with particles inside fl. Since the considered 
system is linear, there exists a Green's function G such that 
Uia = Giajpujp for i £ fl. Let: Ly Q j/3 = Gi a jp — Gj a jp. 
It follows that: Uij a = Lijajpujp. Under a rigid translation 



(all { uj} equal): = 0, hence: ', 
It follows that 



Lijajp [ujp - W/a(r)]. 



Uij a = Lijajp [up(rj) - up(r)] + L i]aJS3 [ujp - up{rj)] , 

(18) 

where ujp — up(vj) is a fluctuating displacement. The sum 
over J in the second term can be shown to be subdominant 
when W sufficiently exceeds w. The first term equals, to lead- 
ing order in a gradient expansion: 



up{rj) -itaO 



(rj 7 - r 7 ) 



(19) 



Substituting Eq. (|l9|) in Eq. (|T|): 

KijLij-yjfi {r Jv — r„) f^j a rijpfij^ (20) 

. ij 

x f dscj>[r - r? + sr%] 
Jo 



o- Q /3(r) 



<>(r), 



where rotational symmetry has been invoked. Thus linear elas- 
ticity is valid when ||e|| <C 1 (the strain components are small) 
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and | WV Q V,gu -c 1. Note that the elastic moduli depend, in 
principle, on the position as well as the resolution (through the 
coarse graining function <f>). Our numerical results (see below) 
indicate that the contribution of the second term on the right 
hand side of Eq. ( |l8[ ) to the stress is smaller than a naive bound 
based on surface to volume ratios implies. 



2.4 Elastic Energy 

In the quasistatic limit, the energy density reduces, cf. Eq. (||), 

e(r,t) = l *(*«(*))#•-*(*)]■ (21) 



to: 



To lowest nonvanishing order in the strain, the potential energy 
corresponding to Eq. (|l5|) is given by \Kij (r°- • uy) .Hence, 
at this order: 

e ( r > ( ) = iE K a (*« ' u ^) 2 ^ ~ >■*(*)]■ ( 22 ) 

If linear elasticity is to hold, one must have: e = i<r ■ e. A set 
of straightforward transformations on Eq. (E2h yields: 



, , 1 d 
e(r '* ) = 2 CT - e -^ 



(23) 

x / dscf)[r - Ti(t) + srij(t)] 
Jo 



It can be shown that the second term on the right hand side of 
Eq. ( |23| ) represents the adiabatic limit of the divergence of the 
heat flux [Eq. (§h], i.e., the work of the fluctuating forces at the 
surface of a control volume. As this term is a divergence of a 
flux, its average over a sufficiently large volume tends to zero, 
i.e. the average of \cr-e over a sufficiently large volume (not its 
'local value') is the elastic energy density. As in the previous 
section, one obtains that classical elasticity is valid only for 
sufficiently large scales, in particular scales for which "surface 
contributions", as explained above, vanish. 



3 Numerical Results 

The above results are demonstrated on a two dimensional (2D) 
system of particles with harmonic interactions [Eq. ([To])]. As a 
first test case, consider a square-shaped triangular lattice con- 
figuration, with uniform nearest-neighbor spring constants K 
and rest lengths equal to the lattice constant d, subjected to the 
following boundary conditions. The displacements of the parti- 
cles at the boundary (whose positions are denoted by { } ) are 
chosen to yield a homogeneous deformation, i.e., the boundary 
is subject to an "applied strain", e a p: 



Ula 



£ a f3 (rip - r p ) , 



(24) 



(by matrix inversion) for a given applied e a p, yielding (for each 
choice of e Q| g) a set of displacements {uf }. The latter are used 
for calculating the linear strain field [using Eq. (14| and linear 
stress, Eq. (|17|). The coarse-graining function used is, 



(R) = <t>(R) = 



A 



R — w/2 

1 + e~ s — 



(25) 



where r^p is an arbitrary point, chosen to reside at a corner of 
the system. The (linear) static equilibrium equations are solved 



(the Fermi distribution), chosen as a smoothed Heaviside func- 
tion. The length on which the function decays to zero, A, can be 
chosen independently of the coarse graining scale, w. In all the 
calculations presented below, A = | . The constant A is fixed 
by the normalization: in 2D, 2n J Q R<f>(R)dR = 1, 

In this case of an ordered lattice, the strain components con- 
verge to the applied strain components even for w = d. This 
should be expected, since, as mentioned, for a lattice configu- 
ration (with uniform spring constants) under homogeneous de- 
formation, the particle displacements correspond to an affine 
transformation, rendering the mean field approximation exact. 
The corresponding stress components are scale-independent as 
well, and consistent with the continuum isotropic elastic mod- 
uli for a triangular lattice (Lame constants A = \i = ^ K ). 

For disordered systems, even for a "homogeneous applied 
strain" as described above, the stress and strain fields are inho- 
mogeneous, in general. Linear elasticity should still be valid (at 
least as an approximation, as described in Sec. |J), i.e. the local 
stress should depend on the local strain by an appropriate lo- 
cal linear relation on a sufficiently large coarse graining scale. 
In order to examine the influence of disorder on the validity 
of linear elasticity, disordered systems were generated, based 
on the triangular lattice configuration: a random number, uni- 
formly distributed in the range [—5d, Sd], is added to the x and 
y coordinates of the particles. Particles whose distance is less 
than c max = l.ld are connected by springs, with rest lengths 
equal to the particle separation in the initial configuration (en- 
suring an unstressed configuration), and spring constants uni- 
formly distributed in the range [K — 6K, K + 8K\. Note that 
this choice of c max can decrease the coordination number, for 
some of the particles, to a value smaller than 6, for sufficiently 
large Sd (this kind of topological disorder can give rise to qual- 
itatively different effects than positional disorder). To evaluate 
the extent to which disordered systems are described by the 
equations of linear elasticity, the following procedure is used. 
Three independent global strains eVg; i = 1,2,3 are applied, 
and the stress and strain fields are calculated at a given point. 
According to linear elasticity, these fields should be linearly 
related, though the elastic moduli may be position dependent. 
In a 2D system with central forces (as used here), there are, in 
general, 6 independent elastic moduli. Each deformation pro- 
vides three linear equations for these moduli (however, it can be 
shown that two independent deformations are insufficient for 
determining them, and three deformations are required). The 
elastic moduli are determined using 6 of the 9 linear equations. 
The stress components which are not used in this procedure 
are then calculated using these elastic moduli, and their values 
are compared to those computed directly using Eq. @. The 
root mean square of the differences between the stress compo- 
nents calculated by employing the measured moduli and their 
directly measured exact values (normalized by the norm of the 
exact values), A, is used as a measure of the extent to which 
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are shown in 



the system is described by linear elasticity at a given position 
and for a given value of the coarse graining scale. 

The stress and strain components at the center of a disor- 
dered system of size 40d x 40d, with Sd = O.ld, 6K = O.IK 

0.005 0.0075 \ 
0.0075 0.01 )' 
Figs. |2||3|, respectively, for different coarse graining scales. These 
fields are obviously inhomogeneous, hence the observed scale 
dependence (note that for coarse graining scales approaching 
the system size, the values of the strain components do ap- 
proach their imposed global values, as expected). 



and applied strain e 1 



viation is quite large for small coarse graining scales, but it 
decreases to less than 1% for w > I2d, indicating that linear 
elasticity holds reasonably well beyond this scale. Note that 
at this scale, the stress and strain are still scale dependent, or 
inhomogeneous, which implies that linear elasticity does hold 
locally inside the system, with elastic moduli which depend on 
position (and scale). The deviations from elasticity for a similar 
system with higher disorder (with 5d — 0.15d, SK = 0.15-fQ 
is shown in Fig. One may interpret the observed fluctuations 
as an indication that the scale required to obtain linear elastic- 
ity exceeds the size of the system, as even at scales close to this 
size there is no clear saturation to a linear elastic relation. 




15 

w/d 

Fig. 2. The strain components vs. the coarse graining scale, w, in a 

disordered configuration (see text), for an applied strain e 1 (see text). Fi §- 4 - The deviation from elasticity in a disordered configuration (see 

text), measured by A (see text), vs. the coarse graining scale, w. 
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Fig. 3. The stress components vs. the coarse graining scale, w, in a 
disordered configuration (see text), for an applied strain e 1 (see text). 



Fig. [I shows the deviations from elasticity, measured by 
A as defined above, for the same disordered system. The de- 
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Fig. 5. The deviation from elasticity in a configuration with higher 
disorder compared to Fig. ^ (see text), measured by A (see text), vs. 
the coarse graining scale, w. 
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4 Concluding Remarks 

It is important to note that the results presented in Sec. |^ have 
been obtained with boundary conditions chosen to obtain a ho- 
mogeneous strain field. Indeed, in the case of a lattice, a uni- 
form strain is obtained. Thus the inhomogeneity in the disor- 
dered case is not a result of the applied boundary conditions 
but of the disorder itself. When inhomogeneous boundary con- 
ditions are applied one expects an inhomogeneous strain field 
even in a homogeneous system. In this case one should observe 
deviations from linear elasticity on small scales or close to the 
boundary. The deviations from standard linear elasticity should 
be particularly prominent in small systems. This is indeed the 
case [|[ for an ordered (lattice configuration) slab of particles 
resting on a rigid support ("bottom") with a point force applied 
to the center particle at the "top" (motivated by experiments 
on granular systems JIo|,pT|]). A comparison of the stress field 
obtained at the bottom of slabs of different heights (number of 
layers of particles) with corresponding continuum elastic so- 
lutions |J shows significant deviations, for a small number of 
layers. These deviations decrease as the size of the system in- 
creases, rendering continuum elasticity a good approximation 
for a sufficient number of layers [O(40) in 2D, O(60) in 3D]. 
For disordered systems, this effect is even more pronounced, as 
elasticity sets in on larger coarse graining scales (which should 
be compared to the size of the system). An additional factor 
which is expected to influence the crossover to linear elastic- 
ity is a possible inhomogeneous stress in the reference state, 
as observed in [|j] (in the examples presented here, the refer- 
ence state is unstressed). In [Q], a similar crossover has been 
observed for the vibrational modes of disordered systems, in 
which the strain is typically inhomogeneous. While not men- 
tioned in fl], their results appear to suggest that the crossover is 
obtained for larger system sizes as the frequency (wavelength) 
is increased, i.e., as the strain gradients are larger, consistently 
with the above arguments. 

The above considerations, as well as the theoretical calcu- 
lations described in Sec. |^, suggest that the general constitu- 
tive description of the systems considered here, even for small 
deformations, should involve a non-local relation between the 
stress and the strain field, a leading approximation thereof be- 
ing provided by gradient elasticity. Larger deformations should 
obviously require nonlinear elasticity. 

While in granular matter, the interactions among the grains 
are not harmonic, we believe that the above discussion may 
still be relevant for quasi-static deformation of granular mate- 
rials. First, the interactions are often described by elastic con- 
tact models for which the force-displacement law, even when 
nonlinear, may be linearized around a reference state (though 
to conform with the nature of cohesionless grains, the springs 
should be "one-sided", i.e., allow for compressive forces only). 
A deformation under which the contact network does not change 
should then still be described by linear elasticity on sufficiently 
large scales. As the interactions among the grains are only com- 
pressive, some contacts may break under a given applied bound- 
ary conditions, yielding a modified contact network. For small 
deformations, the changes in this network may be sufficiently 
small for the elastic moduli not to be affected. When this is 
not the case, incremental elasticity, whereby the elastic mod- 
uli are modified by the deformation, may still be appropriate 



for describing the deformation. Furthermore, when the bound- 
ary conditions result in tensile stress components in a given 
region, contacts may break there preferentially in specific di- 
rections [§[T^]. A similar type of stress-induced anisotropy has 
been suggested in the context of plastic models for soil me- 
chanics [O]. Another source of deviations from elasticity is 
granular friction. Note, however, that static friction is not dis- 
sipative, and it may actually prevent the breaking of contacts, 
extending the elastic range of frictional systems with respect 
to the (idealized) frictionless case. Once the limit of static fric- 
tion is exceeded, friction is kinetic and dissipation does occur. 
In this case one expects plastic failure, which is clearly beyond 
the limits of validity of elastic theory. 

As a final remark, force chains, which have been observed 
in both experiments [14| and simulations [[L5| of granular ma- 
terials, are also observed in elastic systems [Q, ?J. However, as 
shown in [fj], force chains appear even in inhomogeneously 
strained lattices which are macroscopically isotropic, due to the 
inherent small-scale anisotropy of discrete systems. The cor- 
responding stress field, even at small scales, does not exhibit 
similar structures in an isotropic system, i.e., force chains do 
not necessarily imply an inhomogeneous stress field. 

Support from the Israel Science Foundation, grants no. 39/98 and 
53/01, is gratefully acknowledged. 
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